NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


THESIS 

PARAMETRICS  OF  SUBMARINE  DYNAMIC 
STABILITY  IN  THE  VERTICAL  PLANE 

by 

Stavros  1.  Papanikolaou 
March,  1996 

Thesis  Advisor: _ Fotis  A.  Papoulias 

Approved  for  public  release;  distribution  is  unlimited. 


19960530  032 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  0MB  No.  0704 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction,  searching  existing  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other 
ispect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  headquarters  Services,  Directorate  for  Information  Operations  and 
Elcports.  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-01 88) 
Washington  DC  20503. 


1 .  AGENCY  USE  ONLY  (Uave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March  1996  Master's  Thesis 

4.  TITLE  AND  SUBTITLE  Parametrics  of  Submarine  Dynamic  Stability  5.  FUNDING  NUMBERS 
in  the  Vertical  Plane. 


AUTHOR(S):  Stavros  I.  Papanikolaou 


PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 
Naval  Postgraduate  School  REPORT  NUMBER 

Monterey  CA  93943>S000 _ 

SPONSORING/MONITORING  AGENCY  NAME(S)  AND  10.  SPONSORING/MONITORING 

ADDRESS(ES)  AGENCY  REPORT  NUMBER 


1 1 .  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy 
or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT  Approved  for  12b.  DISTRIBUTION  CODE 
public  release;  distribution  unlimited 

13.  ABSTRACT  (maximum  200  words)  The  problem  of  dynamic  stability  of  submersible  vehicles  in  the  dive  plane  is  examined 
utilizing  bifurcation  techniques.  The  primary  mechanism  of  loss  of  stability  is  identified  in  the  form  of  generic  Hopf 
bifurcations  to  periodic  solutions.  Stability  of  the  resulting  limit  cycles  is  established  using  center  manifold  approximations 
and  integral  averaging.  The  hydrodynamic  coefficients  are  calculated  using  existing  semi-empirical  methods.  Parametric 
studies  are  performed  with  varying  vehicle  geometric  properties.  The  methods  described  in  this  work  could  suggest  ways  to 
enlarge  the  submerged  operational  envelope  of  a  vehicle  early  in  the  design  phase. 


14,  SUBJECT  TERMS 

15.  NUMBER  OF  PAGES 

Submarine  stability.  Bifurcations,  Periodic  solutions 

84 

16.  PRICE  CODE 

17.  SECURITY 

CLASSIFICATION  OF 
REPORT 

18.  SECURITY 

CLASSIHCATION 
OF  THIS  PAGE 

19.  SECURITY 

CLASSIHCATION 
OF  ABSTRACT 

20.  LIMITATION  OF  ABSTRACT 
UL 

Unclassified 

Unclassified 

Unclassified 

NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Sid.  239-18 


1 


11 


Approved  for  public  release;  distribution  is  unlimited. 


PARAMETRICS  OF  SUBMARINE  DYNAMIC  STABILITY  IN  THE  VERTICAL 

PLANE 


Stavros  1.  Papanikolaou 
Lieutenant  Jounior  Grade,  Hellenic  Navy 
B.S.,  Hellenic  Naval  Academy,  1989 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 


IV 


ABSTRACT 


The  problem  of  dynamic  stability  of  submersible  vehicles  in  the  dive  plane  is  examined 
utilizing  bifurcation  techniques.  The  primary  mechanism  of  loss  of  stability  is  identified  in  the 
form  of  generic  Hopf  bifurcations  to  periodic  solutions.  Stability  of  the  resulting  limit  cycles  is 
established  using  center  manifold  approximations  and  integral  averaging.  The  hydrodynamic 
coefficients  are  calculated  using  existing  semi-empirical  methods.  Parametric  studies  are 
performed  with  varying  vehicle  geometric  properties.  The  methods  described  in  this  work  could 
suggest  ways  to  enlarge  the  submerged  operational  envelope  of  a  vehicle  early  in  the  design 
phase. 
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1.  INTRODUCTION 


A.  PROBLEM  OVERVIEW 


The  increasing  demands  of  using  submersible  vehicles  for  more  complex  and 
demanding  missions,  force  us  to  use  a  variety  of  methods,  mathematical  mod¬ 
els,  and  assumptions  for  the  study  of  their  dynamic  interactions  and  responses. 
This  study  is  important  in  order  to  enhance  vehicle  operations.  Typically,  lin¬ 
earization  of  the  equations  of  motion  around  nominal  straight  line  level  flight 
paths  along  with  eigenvalue  analysis  can  be  employed  {Arentzen  and  Mandel, 
1960),  (Clayton  and  Bishop,  1982),  (Feldman, 1987).  A  simple  but  efflcient 
stability  criterion  Gv  >  0  can  be  obtained  where  the  stability  index  is 
function  of  the  hydrodynamic  coefficients  in  heave  and  pitch.  Values  for  the 
stability  index  can  be  computed  by. 


M.uj{Zq  -I-  m) 


(1) 


This  index  is  analogous  to  the  familiar  stability  coefficient  for  horizontal  plane 
maneuvering  and  can  be  thought  of  as  a  high  speed  approximation  where  the 
effect  of  the  metacentric  restoring  moment  is  minimal  (Papadimitriou,1994). 
If  the  value  of  Gy  is  greater  than  zero,  the  vehicle  is  dynamically  stable.  As  it 
has  been  established  in  previous  studies  (Papoulias  and  Papadimitriou,  1995) 
though,  this  is  only  a  sufficient,  and  rather  conservative  condition  for  stability. 
Nevertheless,  it  is  widely  used  and  its  value  is  indicative  of  vertical  plane  sta¬ 
bility  for  any  new  design.  We  should  keep  in  mind,  however,  that  the  condition 
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(?v  <  0  indicates  a  divergent  loss  of  stability  which  is  quite  uncommon  in  the 
vertical  plane.  Most  modern  submarines  exhibit  a  flutter-like  instability  at 
high  speed,  which  can  not  be  analyzed  using  the  above  simplified  index.  Di¬ 
vergent  motions  may  develop  in  combined  six  degrees  of  freedom  (Papoulias  et 
al,  1993)  and  their  occurrence  can  not  be  analyzed  by  a  single  stability  index. 
Previous  work  (Papadimitriou,  1994)  was  limited  to  a  single  body  with  fixed 
hydrodynamic  coefficients.  In  this  work,  we  expand  by  allowing  the  geometry 
of  the  body  and  thus  its  hydrodynamic  properties  to  vary. 


B.  THESIS  OUTLINE 

Previous  work  (Papoulias  and  Papadimitriou,  1995)  analyzed  the  problem 
of  stability  of  motion  with  controls  fixed  in  the  vertical  plane,  with  partic¬ 
ular  emphasis  on  the  mechanism  of  loss  of  stability  of  straight  line  motion. 
The  closed  loop  control  problem  was  analyzed  in  (Papoulias  et  al,  1995).  The 
surge  equation  was  decoupled  from  heave/pitch  through  a  perturbation  series 
approach  (Bender  and  Orszag,  1978).  As  was  established  in  (Papadimitriou, 
1994)  loss  of  stability  occurs  in  the  form  of  generic  bifurcations  to  periodic 
solutions  (Guckenheimer  and  Holmes,  1983).  Taylor  expansions  and  center 
manifold  approximations  were  employed  in  order  to  isolate  the  main  nonlinear 
terms  that  influence  system  response  after  the  initial  loss  of  stability  (Hassard 
and  Wan,  1978).  Integral  averaging  was  performed  in  order  to  combine  the 
nonlinear  terms  into  a  design  stability  coefficient  (Chow  and  Mallet-Paret, 
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1977) .  Some  difficulties  associated  with  the  nonsmoothness  of  the  absolute 
value  nonlinearities  was  dealt  with  by  employing  the  concept  of  generalized 
gradient  (Clarke,  1983).  This  was  employed  as  an  alternative  to  the  lin¬ 
ear/cubic  approximation  typically  used  in  ship  roll  motion  studies  (Dalzell, 

1978) .  The  same  methodology  is  applied  in  this  work  in  order  to  analyze  the 
sensitivity  of  the  results  with  respect  to  geometric  characteristics  of  the  body. 

Vehicle  modeling  in  this  work  follows  standard  notation  (Gertler  and  Ha¬ 
gen,  1976),  (Smith  et  al,  1978),  and  numerical  results  are  presented  for  a  family 
of  bodies  of  revolution  similar  to  the  DARPA  SUBOFF  model  (Roddy,  1990) 
for  which  a  set  of  hydrodynamic  coefficients  and  geometric  properties  is  avail¬ 
able.  This  parametric  study  is  conducted  utilizing  existing  senai-empirical 
methods  for  the  calculation  of  hydrodynamic  coefficients.  The  methods  are 
based  on  (Fidler  and  Smith,  1978),  (Humphreys  and  Watkinson,  1978),  (Pe¬ 
terson,  1980)  and  have  been  verified  in  (Wolkerstorfer,  1995).  The  effects  of 
varying  the  nose,  base,  and  tail  fractions  of  the  body  as  well  its  nondimen- 
sional  volume  to  length  ratio  on  the  hydrodynamic  derivatives  were  studied  in 
(Holmes,  1995)  where  prediction  equations  were  derived  based  on  curve  fitting 
of  the  results.  These  hydrodynamic  prediction  equations  are  normalized  by 
taking  the  SUBOFF  model  as  a  baseline.  This  model  has  been  experimentally 
validated  for  angles  of  attack  on  the  hull  between  ±15  deg.,  while  the  constant 
coefficient  approximation  introduces  very  little  error  in  time  domain  simula¬ 
tions  (Tinker,  1978).  Unless  otherwise  mentioned,  all  results  in  this  work  are 
presented  in  standard  dimensionless  form  with  respect  to  the  vehicle  length 
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L  =  4.26  m,  and  nominal  forward  speed  U  =  2.44  m/sec  (Papadimitrion, 
1994). 
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II.  PROBLEM  FORMULATION 


A.  EQUATIONS  OF  MOTION 

In  order  to  obtain  the  mathematical  model  the  following  assumptions,  re¬ 
strictions,  and  definitions  have  to  be  made: 

1.  The  submersible  vehicle  motion  is  restricted  in  the  vertical  plane,  thus 
the  model  consists  of  coupled  nonlinear  heave  and  pitch  equations. 

2.  The  coordinate  frame  is  fixed  at  the  vehicle’s  geometrical  center. 

3.  Vehicle  is  port/starboard  symmetric  and  neutrally  buoyant. 

4.  Use  Newton’s  equations  of  motion  in  dimensionless  form. 

The  nonlinear  heave  and  pitch  equations  become: 

m{w  —  uq  —  zaq^  -  xcq)  =  Zqq  Z^w  +  Zqq  ZyjW 

rnose 

—Cd  /  h{x){w  —  xq)\w  —  xq\dx,  (2) 

Jtail 

lyq  -I-  mza{u  +  wq)  —  ruxc^w  —  uq)  =  Mqj  4-  M^w  +  Mqq  +  MyjW 

rnose 

+Cz)  /  b{x){w  —  xq)\w  —  xq\x  dx 

Jtzil 

— xgbW  cos  0  —  zqbW  sin  9,  (3) 

where  xqb  =  xq  —  xb,  zgb  =  zq  —  zb,  and  the  rest  of  the  symbols  are  based 
on  standard  notation  as  shown  in  Table  1.  Without  loss  of  generality  we  can 
assume  that  zb  =  2:5  =  0,  so  that  xqb  —  ^gb  =  zq-  The  cross  flow 
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integral  terms  in  these  equations  become  very  important  for  high  angles  of 
attack  maneuvering,  where  they  provide  the  primary  motion  damping.  The 
drag  coefficient,  Co,  is  assumed  to  be  constant  throughout  the  vehicle  length 
for  simplicity.  This  does  not  affect  the  qualitative  properties  of  the  results 
that  follow.  The  vehicle  pitch  rate  is, 

»  =  ,.  (4) 

Dynamic  coupling  between  surge  and  heave /pitch  is  present  due  to  coordinate 
coupling  as  a  result  of  the  nonzero  metacentric  height.  However,  it  has  been 
shown  (Papoulias  and  Papadimitriou,  1995)  that  this  coupling  is  of  higher 
order  and  does  not  change  the  linear  and  nonlinear  results  that  follow. 


B.  HYDRODYNAMIC  COEFFICIENTS 

Systematic  studies  based  on  semi-empirical  methods  have  resulted  in  the 
evaluation  of  hydrodynamic  coefficients  for  a  generic  body  of  revolution  in 
terms  of  basic  geometric  properties.  Curve  fitting  revealed  that  adequate  ac¬ 
curacy  for  initial  design  can  be  obtained  by  equations  of  the  form 

He  =  AiF^  +  A2FnFm  +  AsFl  +  A4Fn 

d-AsEm  +  Ae  +  Ay  , 

where  He  denotes  a  given  coefficient  in  its  standard  nondimensional  form,  V 
the  underwater  volume  of  the  body,  L  its  nominal  length,  the  nose  fraction, 
and  Fm  the  mid-body  fraction.  The  regression  coefficients  A,  are  presented 
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Ai 

regression  coefficient 

b(x) 

local  beam  of  the  hull 

C 

nominal  value  of  volumetric  coefficient 

Cd 

quadratic  drag  coefficient 

Fn 

nose  length  fraction 

middle-body  length  fraction 

He 

given  hydrodynamic  coefficient 

4 

vehicle  mass  moment  of  inertia 

K 

nonlinear  stability  coefficient 

L 

vehicle  length 

m 

vehicle  mass 

M 

pitch  moment 

Ma 

derivative  of  M  with  respect  to  a 

Q 

pitch  rate 

T 

transformation  matrix  of  x  to  z 

U 

forward  speed 

Uq 

critical  value  of  u 

V 

total  volume 

w 

heave  velocity 

X 

state  variables  vector,  x  =  [6,w,q] 

{xb,  Zb) 

body  fixed  coordinates  of  vehicle  center  of  buoyancy 

{XG,  zg) 

body  fixed  coordinates  of  vehicle  center  of  gravity 

z:gb 

center  of  gravity /center  of  buoyancy  separation,  xg  —  xb 

zgb 

vehicle  metacentric  height,  zq  —  zb 

z 

heave  force 

Za 

derivative  of  Z  with  respect  to  a 

OCij 

expansion  coefficients  of  zs  in  terms  of  zi,  z<i 

6 

stern  plane  deflection 

e 

criticality  difference,  e  —  u  —  Uc 

e 

pitch  angle 

Table  1:  Nomenclature 
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He 

Ai 

^2 

^3 

A4 

^6 

Aj 

Zyj 

-0.0641 

-0.1149 

-0.0632 

+0.0670 

+0.0732 

-0.0263 

-0.5769 

+0.0277 

+0.0499 

+0.0266 

-0.0283 

-0.0301 

Zq 

-0.0314 

-0.0559 

-0.0292 

-0.0091 

-0.0880 

■n 

-0.0003 

+0.0040 

+0.0027 

+0.0006 

-0.1590 

+0.0002 

+0.0007 

+0.0007 

-0.0008 

-0.0016 

-0.0144 

-1.8067 

-0.0002 

-0.0007 

-0.0007 

+0.0144 

+1.8067 

wm 

+0.0031 

+0.0024 

-0.0013 

-0.0808 

Table  2:  Regression  coefficients  Ai 


in  Table  2.  Zq  was  assumed  constant  since  the  semi-empirical  techniques 
failed  to  compute  a  reliable  value.  Basic  geometric  definitions  for  the  body 
are  presented  in  Figure  1.  The  constant  C  is  approximately  8  x  10“^  and  is 
the  nominal  value  for  the  volumetric  coefficient.  These  expressions  are  for  a 
body  of  revolution  without  appendages  and  assume  parabolic  nose,  parallel 
mid-body,  and  conical  tail  (Holmes,  1995).  Typical  ranges  of  applicability  for 
these  regression  formulas  are  0.05  to  0.25  for  Fn,  0.40  to  0.60  for  Fm,  and  6.0  to 
10.0  for  V/L^.  Sample  results  for  the  above  hydrodynamic  coefficients  versus 
the  nose  and  mid-body  fraction  ratios  are  presented  in  Figures  2  through  8. 


C.  DEGREE  OF  STABILITY 

The  degree  of  stability  is  defined  as  the  largest  real  part  of  all  eigenvalues  of 
the  linearized  system  of  equations  (2),  (3),  and  (4).  Positive  values  indicate  an 
unstable  system  while  negative  values  show  stability  of  forward  motion.  The 
degree  of  stability  versus  xqb  for  constant  forward  speed  u  =  0.5  and  different 
values  of  zgb  is  shown  in  Figures  9  through  12.  Based  on  these  results  we  can 
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zwdot 


xIO"* 


Figure  3:  Hydrodynamic  coefficient  Mw  versus  Fn  and  Fm 


Figure  4:  Hydrodynamic  coefficient  versus  Fn  and  Fm 
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12 


Figure  9:  Degree  of  stability  for  u  =  0.5,  varying  zqb^  Fn  =  0.3,  and  Pm  —  0'6 
draw  the  following  conclusions: 

1.  In  all  cases  the  vehicle  is  dynamically  more  stable  as  the  metacentric 
height  zqb  is  increased. 

2.  In  all  cases  the  vehicle  is  dynamically  less  stable  as  the  separation  between 
the  centers  of  gravity  and  buoyancy  is  reduced  in  absolute  value. 

3.  For  constant  Fn,  increasing  values  of  Fm  result  in  less  stable  vehicles. 
This  means  that  a  longer  tail  is  beneficial  for  stability  of  motion,  as 
expected. 

4.  The  same  conclusion  holds  for  constant  mid-body  ratio  and  varying 
nose  ratios  Fn- 
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Figure  11:  Degree  of  stability  for  u  =  0.5,  varying  zqb^  Fn  —  0.3,  and 


Figure  12:  Degree  of  stability  for  u  =  0.5,  varying  zgb,  Fn  =  0.1,  and  Fm  =  0.6 

Corresponding  results  for  constant  zqb  =  0.015  and  varying  forward  speeds 
u  are  shown  in  Figures  13  through  16.  Similar  conclusions  as  those  discussed 
previously  hold  in  these  cases  with  the  following  exceptions: 

1.  For  very  low  forward  speeds,  the  case  xg  =  0  may  be  best  for  stability. 

2.  For  very  low  speeds,  smaller  tails  may  result  in  more  stable  configura¬ 
tions. 

Combined  results  for  variations  in  both  xqb  and  u  are  shown  by  the  mesh 
plots  of  Figures  17  through  20.  The  value  of  zqb  was  held  constant  at  0.015 
for  all  plots.  These  figures  confirm  our  previous  conclusions  by  presenting  the 
results  in  more  detail. 

Figure  21  shows  the  degree  of  stability  versus  Fn  and  Fm-  Both  values 
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Figure  17:  Degree  of  stability  for  Fn  =  0.3  and  Fm  =  0.6 
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Figure  19:  Degree  of  stability  for  Fn  =  0.3  and  Fm  —  0.4 
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of  iCG  and  zq  were  kept  constant  and  equal  to  0  and  0.015  respectively.  The 
three  surfaces  shown  correspond  to  values  u  =  0.4, 0.5, 0.6.  The  upper  one 
corresponds  to  u  =  0.6  while  the  lower  one  to  u  =  0.4.  It  can  be  seen  that 
the  degree  of  stability  becomes  more  negative  for  decreasing  u,  and,  generally 
speaking,  for  decreasing  Fn  and  Fm- 

Figure  22  shows  the  degree  of  stability  versus  Fn  and  Fm-  Both  values  of 
forward  speed  u  and  zq  were  kept  constant  and  equal  to  0.5  and  0.015  respec¬ 
tively.  The  three  surfaces  shown  correspond  to  values  xg  =  -0.01,0,  +0.01. 
The  upper  one  corresponds  to  xg  =  0.0  while  the  lower  one  to  xg  =  +0.01.  It 
can  be  seen  that  the  degree  of  stability  becomes  more  negative  for  increasing 
xg  in  absolute  value,  and,  generally  speaking,  for  decreasing  Fn  and  Fm- 

Figure  23  shows  the  degree  of  stability  versus  Fn  and  Fm-  Both  values  of  for¬ 
ward  speed  u  and  xg  were  kept  constant  and  equal  to  0.5  and  0.0  respectively. 
The  three  surfaces  shown  correspond  to  values  zg  =  +0.005, +0.015, +0.025. 
The  upper  one  corresponds  to  zg  =  +0.005  while  the  lower  one  to  zg  = 
+0.0025.  It  can  be  seen  that  the  degree  of  stability  becomes  more  negative  for 
increasing  zg,  and,  generally  speaking,  for  decreasing  Fn  and  Fm- 


D.  CRITICAL  SPEED 

The  parameter  value  where  the  real  part  of  the  dominant  complex  conjugate 
pair  of  eigenvalues  crosses  zero  defines  the  point  where  linear  stability  is  lost. 
This  critical  point  can  be  computed  by  considering  the  characteristic  equation 
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Figure  21:  Degree  of  stability  versus  Fn  and  Fm  for  xq  =  0,  zg  =  0.015,  and 
u  =  0.4, 0.5, 0.6 


0.35 


Figure  22:  Degree  of  stability  versus  F„  and  Fm  for  u  =  0.5,  zq  =  0.015,  and 

XQ  =  —0.01, 0,  +0.01 


21 


Figure  23:  Degree  of  stability  versus  Fn  and  Fm  for  u  =  0.5,  xq  =  0,  and  zg  - 
0.005,0.015,0.025 


of  the  system  (Papadimitriou,  1994).  Routh’s  criterion  applied  to  this  can  be 
solved  for  the  dimensionless  weight. 


W  = 


^2^2,0 

A2D2,\  —  B2C2,1 


(5) 


where. 


C2,0  “  Zxj}{Mq  —  mxg)  —  MiuiZq  +  Tn)  , 

C2,\  —  {m-  Zw)  {zGB  cos  9q  -  XQB  sin  0o)  , 
£>2,1  =  .^^(a^GBsin^o  -  ■^gbcos^o)  • 


It  should  be  mentioned  that  the  effect  of  the  forward  speed  u  is  embedded  into 
the  definition  for  the  dimensionless  vehicle  weight  W  through, 


(6) 
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The  value  of  the  critical  speed  Uc  can  then  be  evaluated  from  (5)  and  (6). 
Typical  results  are  presented  in  Figures  24  through  27.  A  family  of  critical 
speeds,  «c,  is  shown  versus  xq  with  zq  as  the  parameter  of  the  curves.  These 
results  were  obtained  for  a  nose  fraction  F„  =  0.1, 0.3  and  mid-body  fraction 
Fm  =  0.4, 0.6.  The  volumetric  coefficient  was  kept  at  nominal  for  all  results. 
Vertical  plane  motions  are  stable  for  forward  speeds  less  than  the  critical  speed. 
It  can  be  seen  that  stability  is  increasing  with  increasing  zq  while  xg  =  0  is  the 
most  conservative  condition  for  stability.  Therefore,  a  vehicle  which  is  stable 
when  properly  trimmed  will  remain  stable  for  off-trim  conditions.  The  fact 
that  a  vehicle  with  a  longer  aft-body  ought  to  be  dynamically  more  stable  is 
confirmed  by  comparing  the  results  of  Figures  24  and  26  to  the  results  shown 
in  Figures  25  and  27  respectively.  It  can  be  seen  that  the  corresponding  critical 
speeds  become  smaller,  thereby  reducing  the  dynamic  stability  margin,  as  the 
nose  and  mid-body  fractions  are  raised.  This  trend  is  consistent  for  all  values 
of  xg  and  zq  examined. 

Combined  plots  of  the  critical  speed  versus  both  xq  and  zq  are  shown  in 
Figures  28  and  29.  Figure  28  presents  the  surfaces  for  Fn  =  0.3  and  Fm  = 
0.4, 0.5, 0.6.  The  uppper  surface  corresponds  to  Fm  =  0.4.  Figure  29  presents 
the  surfaces  for  Fm  =  0.5  and  Fn  =  0.1, 0.2, 0.3.  The  upper  surface  corresponds 
to  Fn  =  0.1. 

Combined  plots  of  the  critical  speed  versus  both  Fn  and  Fm  are  shown 
in  Figures  30  through  32.  Figure  30  presents  the  surface  when  zq  =  0.0125 
and  Xq  =  0.  Figure  31  gives  us  a  comparative  view  keeping  zg  =  0.0125  and 
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Figure  26:  Critical  speed  versus  xq  for  Fn  =  0.3  and  Fm  =  0.4  and  different  values 
of  ZG 


Figure  27:  Critical  speed  versus  xq  for  Fn  =  0.3  and  Fm  =  0.6  and  different  values 


Figure  28:  Critical  speed  versus  xg  and  zg  for  Fn  =  0.3  and  Fm  =  0.4, 0.5, 0.6 


Figure  29:  Critical  speed  versus  xq  and  zg  for  Fm  =  0.5  and  Fn  =  0.1, 0.2, 0.3 


26 


28 


using  XG  =  —0.01,0, +0.01  to  plot  the  surfaces  as  shown.  The  lower  surface 
corresponds  to  xq  =  0.  It  can  be  seen  that  nonzero  xg  increases  the  range 
of  stability,  while  the  general  trend  is  to  increase  stability  as  both  Fn  and 
Fm  become  smaller.  A  similar  plot  for  xg  =  0  and  for  three  values  of  zg, 
0.005,  0.010,  and  0.025  is  shown  in  Figure  32.  The  lower  surface  corresponds 
to  ZG  =  0.005  and  the  higher  one  to  zg  —  0.025.  It  can  be  seen  that  the 
metacentric  height  has  by  far  the  greatest  effect  on  dynamic  stability,  while 
the  effects  of  hull  geometry  are  smaller. 

For  comparison,  a  plot  of  the  classical  stability  coefficient  Gv  from  equation 
(1),  is  shown  in  Figure  33.  The  different  curves  correspond  to  various  mid¬ 
body  fractions,  while  the  nose  fraction  is  kept  constant.  It  can  be  seen  that  Gy 
is  negative  throughout.  Therefore,  it  would  have  predicted  an  unstable  vehicle 
for  all  ranges  of  the  parameters,  which  is  of  course  incorrect.  Furthermore,  Gt, 
becomes  less  negative  as  F^  is  increased,  which  would  suggest  that  dynamic 
stability  is  increased  as  the  aft-body  length  is  decreased.  This  is  also  a  false 
conclusion.  As  we  pointed  out  in  the  introduction,  the  classical  stability  index 
G„  should  be  used  with  extreme  caution. 
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III.  BIFURCATION  ANALYSIS 


A.  INTRODUCTION 

The  nonlinear  bifurcation  analysis  is  based  on  the  general  methodology 
used  in  (Papadimitriou,  1994).  The  fundamental  equations  are  reproduced 
here  for  completeness  of  the  presentation.  The  nonlinear  heave/pitch  equations 
of  motion  (2),  (3),  and  (4)  are  written  in  the  form, 

e  =  q,  (7) 

w  =  aiiw  +  ai2q  +  ai3{xGB  cos  6  +  zqb  sin  0) 

+dy,{w,q)  +  ci{w,q)  ,  (8) 

q  =  a2iw  +  a22q  +  a23(a^C?B  cos  6  +  zqb  sin  6) 

+dq{w,q) +  C2iw,q)  ,  (9) 

where  the  various  coefficients  are  functions  of  the  hydrodynamic  derivatives 
and  mass  properties,  and  7^,,  I,  are  the  cross  flow  integrals. 

The  system  of  equations  (7)  through  (9)  is  written  in  the  compact  form 

x  =  Ax  +  g(x),  (10) 

where 

x=[0,io,g],  (11) 

is  the  three  state  variables  vector,  and  A  is  the  linearized  sytem  matrix  eval¬ 
uated  at  the  nominal  point  xq.  The  term  g(x)  contains  all  nonlinear  terms 
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of  the  equations.  Hopf  bifurcation  analysis  can  be  performed  by  isolating  the 
primary  nonlinear  terms  in  g(x).  Keeping  terms  up  to  third  order,  we  can 
write 

g(x)  =  g^^^(x)  +  g^^^x)  .  (12) 

Using  equations  (7)  through  (11),  the  various  terms  in  (12)  can  be  written  as, 


=  0, 

92^  =  (4  “  M^)mzGq^  -  [mxQ  +  Zq)mzGwq 

9)  >  (13) 

=  —{m  —  Z^)mzGwq  +  {mxa  +  M.w)mzGq^ 
q)  , 

=  0, 

5^2^^  =  d^^\w,q)  + 

Io-iz{xgb  sin  60  -  zgb  cos  ^o)^^  ,  (14) 

93^  =  df\w,q)  + 

lo,23ixGB  sin  00  -  zgb  cos  0o)9^  • 


Expansion  in  Taylor  series  of  dg  requires  expansion  of  the  cross  flow  inte¬ 
grals  1^,  Iq,  which  require  the  Taylor  series  of 

=  (15) 


This  expression  can  be  converted  into  an  analytic  function  using  Dalzell’s 
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approximation  (Dalzell,  (1978), 


48^0  ’ 

which  is  derived  by  a  least  squares  fit  of  an  odd  series  over  some  assumed  range 
of  namely  — <  ^  <  ^c-  This  approximation  has  been  extensively  used  in 
ship  roll  motion  studies  and  is  very  useful  for  its  intended  purpose.  However, 
in  the  present  problem  it  suffers  from  the  several  drawbacks  (Papadimitriou, 
1994).  Instead  of  Dalzell’s  approximation,  we  employ  the  concept  of  general¬ 
ized  gradient  (Clarke,  (1983),  which  is  used  in  the  study  of  control  systems 
involving  discontinuous  or  non-smooth  functions.  In  this  way  we  approximate 
the  gradient  of  a  non-smooth  function  at  a  discontinuity  by  a  map  equal  to  the 
convex  closure  of  the  limiting  gradients  near  the  discontinmty.  In  our  problem 
we  write. 


/(O  =  6l^o|  +  2|^o|(^  -  ^o)  + 

sign(fo)({-&f +  /'’>{«),  (17) 

as  the  Taylor  series  epansion  of  /(^)  near  ^o-  The  sign  function  in  (17)  can  be 
approximated  by, 

sign(^o)  =  lim  tanh  ( —  )  .  (18) 

7->0  \1  J 

The  quantity  7  is  a  small  regularization  parameter  and  is  used  for  proper 
normalization  of  the  results.  Using  (18),  we  can  approximate  /(^)  in  the 
vicinity  of  ^0  =  0  by, 

?i?i «  ■  (19) 
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Since 


^  w  —  xq  ,  (20) 

we  can  express  the  non-smooth  cross  flow  integral  terms  by, 

Iw  =  ^{EqW^  -  SEiW^q  +  3E2wq^  -  E^q^)  , 

67 

Ig  =  ^{Eiw^  -  3E2w\  +  SEswq'^  -  E^q^)  , 

67 

where 


are  the  moments  of  the  vehicle  “waterplane”  area. 

Using  the  previous  second  and  third  order  Taylor  series  expansions,  equa¬ 
tion  (10)  is  written  in  the  form, 

X  =  Ax  -4-  g^^^(x)  -H  g^^^(x)  .  (22) 


If  T  is  the  matrix  of  eigenvectors  of  A  evaluated  at  the  critical  point  u  —  Uc, 
the  linear  change  of  coordinates, 

X  =  Tz  ,  z  =  T"^x  ,  (23) 


transforms  system  (22)  into  its  normal  coordinate  form. 


z  =  T"^ATz  +  T”^g^^)(Tz)  +  T~^g(^)(Tz)  . 

At  the  Hopf  bifurcation  point,  matrix  T“^AT  takes  the  form, 

0  — cjQ  0 

T"^AT  =  I  wo  0  0 

0  0  p 


(24) 
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where  uq  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 
remaining  eigenvalue  p  is  negative.  For  values  of  u  close  to  the  bifurcation 
poit  Uc,  matrix  T“^AT  becomes, 


T"^AT  = 


a'e  —(ujo  +  Lo'e)  0 
(ll>o  +  uj'e)  a'e  0 

0  0  p  +  p'e 


where  e  denotes  the  criticality  difference 


e  =  u  —  Uc , 


and 


(25) 


a'  =  derivative  of  the  real  part  of  the  critical 
eigenvalue  with  respect  to  e  , 
u'  =  derivative  of  the  imaginary  part  of  the 
critical  eigenvalue  with  respect  to  e  , 
p'  =  derivative  of  p  with  respect  to  e  . 

Due  to  continuity,  the  eigevalue  p+p'e  remains  negative  for  small  nonzero 
values  of  e.  Therefore,  the  coordinate  corresponds  to  a  negative  eigenvalue 
and  is  asymptotically  stable.  Center  manifold  theory  predicts  that  the  rela¬ 
tionship  between  the  critical  coordinates  zi,  Z2  and  the  stable  coordinate  zz  is 
at  least  of  quadratic  order.  We  can  then  write  zz  as, 

^3  =  oiuzl  +  anziz2  +  0:22^  >  (26) 

where  the  coefficients,  in  the  quadratic  center  manifold  expansion  (26) 
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need  to  be  determined.  By  differentiating  equation  (26)  we  obtain, 

23  =  2aii2iii  +  Q:i2('^1^2  +  ^1^2)  +  2q:22'22'22  ■  (27) 

We  substitute  ii  =  — a;0'22  ^2  —  we  obtain 

—  0^12^0^!  +  2(q!22  ~  Q!1i)^0'2i^2  ~  Q!l2^0^2  ■  (28) 

The  third  equation  of  (24)  is  written  as, 

i3  =  p23+[T-V^)(Tz)]^^^g^  ,  (29) 

where  terms  up  to  second  order  have  been  kept.  If  we  denote  the  elements  of 
T  and  by, 

T  =  [rriij]  ,  T-^  =  K]  ,  (30) 

then 

X-ig(2)(Tz)  = 

where  expressions  for  di,  ^2,  ds,  and  the  coefficients  iij  are  given  in  Papadim- 
itriou  (1994). 

Equation  (29)  then  becomes 

h  =  pzz  +  ds  ,  (31) 

and  substituting  (26)  and  the  expression  for  dz  into  (31)  we  get, 

23  =  (po^ii  +  ^32^25  +  ^33^35)21 

+  (pai2  +  nz2^26  +  ?^33^36)2i22 
+  (P^22  +  ^32^27  +  ^33^37)22  •  (82) 
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Comparing  coefficients  of  (28)  and  (32)  we  get  a  system  of  linear  equations 
which  yields  the  coefficients  in  the  center  manifold  expansion  (26). 

Using  the  previous  Taylor  expansions  and  center  manifold  approximations, 
we  can  write  the  reduced  two-dimensional  system  that  describes  the  center 
manifold  flow  of  (24)  in  the  form, 

ii  =  a' ezi  -  (wo  +  oj'€)z2  +  ^1(21, 22)  , 
h  =  (t^O  +  oj'e)zi  +  a'€Z2  +  ^2(21, 22)  , 

where  Fi,  F2  are  cubic  polynomials  in  21  and  22. 

If  we  introduce  polar  coordinates  in  the  form, 

zi  =  Rcos<p  ,  22  =  i? sin </» , 

we  can  produce  an  equation  describing  the  rate  of  change  of  the  radial  coor¬ 
dinate  R, 

R  =  a'eR  +  Pi<j))R^  -b  Qi<p)R^  . 

This  equation  contains  one  variable,  R,  which  is  slowly  varying  in  time,  and 
another  variable,  <^,  which  is  a  fast  variable.  Therefore,  it  can  be  averaged 
over  one  complete  cycle  in  (f>  to  produce  an  equation  with  constant  coefficients 
and  similar  stability  properties, 

it  =  a'eR  -b  KR^  +  LR^  , 

where 
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=  g(3rii  +  ri3  +  r22  +  3r24)  , 

1  /■2’f 

L  =  IT  = 

ZTT  70 

Therefore,  the  averaged  equation  becomes 

R  =  a’eR  +  KR^  .  (33) 

Equation  (33)  admits  two  steady  state  solutions,  one  at  i?  =  0  which 
corresponds  to  the  trivial  equilibrium  solution  at  zero,  and  one  at 

Ro  =  .  (34) 

This  equilibrium  solution  corresponds  to  a  periodic  solution  or  limit  cycle  in 
the  cartesian  coordinates  zi,  z^.  For  this  limit  cycle  to  exist,  the  quantity  Rq 
must  be  a  real  number.  In  our  case  a'  is  always  positive,  since  the  system  loses 
its  stability;  i.e.,  the  real  part  of  the  critical  pair  of  eigenvalues  changes  from 
negative  to  positive,  for  increasing  u.  Therefore,  existence  of  these  periodic 
solutions  depends  on  the  value  of  K.  Specifically, 

•  if  iiT  <  0,  periodic  solutions  exist  for  e  >  0  or  u  >  Uc,  and 

•  K  >0,  periodic  solutions  exist  for  e  <  0  or  u  <  Uc- 

The  characteristic  root  of  (33)  in  the  vicinity  of  (34)  is 

j3  =  -2a' e  ,  (35) 


and  we  can  see  that 
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Figure  34:  Nonlinear  stability  coefficient  versus  xq  for  =  0.1,  Fm  =  0.4,  and 
different  values  of  Cd 

•  if  periodic  solutions  exist  for  u>  Uc  they  are  stable,  and 

•  if  periodic  solutions  exist  for  u  <  they  are  unstable. 


B.  RESULTS  AND  DISCUSSION 

Typical  results  of  the  nonlinear  stability  coefficient  K  are  shown  in  Figures 
34  through  37.  Figure  35  presents  a  plot  of  Tf  •  7  versus  xq  for  zq  =  0.015, 
F„  =  0.1,  Fm  =  0.6,  and  for  different  values  of  the  quadratic  drag  coefficient 
Ci).  It  should  be  emphasized  that  the  use  of  JiT  -7  is  more  meaningful  than  the 
use  of  F,  since  it  properly  accounts  for  the  use  of  the  regularization  parameter 
7.  Numerical  evidence  suggests  that  all  curves  K  ■  7  versus  xq  converge  for 
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Figure  35:  Nonlinear  stability  coefficient  versus  xq  for  Fn  =  0.1,  Fm  =  0.6,  different 
values  oi  Cd 


Figure  36:  Nonlinear  stability  coefficient  versus  xq  for  Fn  =  0.3,  Fm  =  0.4,  different 
values  of  C  n 


Figure  37:  Nonlinear  stability  coefficient  versus  xq  for  Fn  =  0.3,  Fm  =  0.6,  different 
values  of  Cp 

7  — >  0.  For  practical  purposes,  values  of  j  smaller  than  0.001  produce  identi¬ 
cal  results.  The  results  of  Figure  8  demonstrate  the  profound  effect  that  the 
quadratic  drag  coefficient  Cp  has  on  stability  of  limit  cycles.  All  Hopf  bifur¬ 
cations  are  supercritical  (K  <  0),  and  they  become  stronger  supercritical  as 
Cp  is  increased.  It  is  worth  noting  that  results  for  Cp  =  0  produce  subcritical 
behavior,  K  >  0,  which  is  clearly  incorrect.  Thus,  neglecting  the  effects  of  Cp 
would  have  produce  entirely  wrong  results  in  the  present  problem.  Additional 
results  show  that  the  bifurcations  become  stronger  supercritical  as  initial  sta¬ 
bility  ZQ  is  increased.  Figure  34  presents  similar  results  with  the  only  difference 
being  the  value  of  mid-body  fraction  Fm  =  0.4.  It  can  be  seen  that  smaller 
Fm  for  the  same  Fn,  which  results  in  longer  body  tail,  may  be  beneficial  for 
stability  in  the  linear  sense  but  it  also  generates  less  supercritical  bifurcations. 
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Figure  38:  Nonlinear  stability  coefficient  versus  Fm  for  =  0,  Cx)  =  0.5,  and 
different  values  of  Fn 

This  can  probably  be  attributed  to  the  increased  responsiveness  of  the  vehi¬ 
cle.  Figures  36  and  37  show  the  same  results  for  nose  fraction  Fn  =  0.3.  It 
should  be  emphasized,  however,  that  altering  the  fore  and  aft  body  lengths 
might  infuence  the  values  of  Cn  which,  as  we  pointed  out,  is  the  single  most 
important  parameter  for  the  nonlinear  nature  of  the  bifurcations. 

Figure  38  shows  the  nonlinear  stability  coefficient  versus  Fm  for  different 
values  of  F„,  while  xq  =  ^  and  Cd  =  0.5.  It  can  be  seen  that  smaller  Fn 
for  the  same  Fm^  which  results  in  longer  body  tail,  generates  less  supercritical 
bifurcations. 

Figure  39  shows  tne  nonlinear  stability  coefficient  versus  Fn  for  different 
values  of  Fm-,  while  xq  =  ^  and  Cd  =  0.5.  Again  it  is  clear  that  longer  body 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  work  presented  a  comprehensive  nonlinear  study  of  straight  line  sta¬ 
bility  of  motion  of  submersibles  in  the  dive  plane  under  open  loop  conditions. 
A  systematic  perturbation  analysis  demonstrated  that  the  effects  of  surge  in 
heave/pitch  are  small  and  can  be  neglected.  Primary  loss  of  stability  was 
shown  to  occur  in  the  form  of  Hopf  bifurcations  to  periodic  solutions.  The 
critical  speed  were  instability  occurs  was  computed  in  terms  of  met  acentric 
height,  longitudinal  separation  of  the  centers  of  buoyancy  and  gravity,  and 
the  dive  plane  angle.  Analysis  of  the  periodic  solutions  that  resulted  from  the 
Hopf  bifurcations  was  accomplished  through  Taylor  expansions,  up  to  third 
order,  of  the  equations  of  motion.  A  consistent  approximation,  utilizing  the 
generalized  gradient,  was  used  to  study  the  non-analytic  quadratic  cross  flow 
integral  drag  terms.  The  main  results  of  this  study  are  summarized  below: 

1.  The  critical  speed  of  loss  of  stability  is  a  monotonically  increasing  func¬ 
tion  of  both  vertical  and  longitudinal  LCG/LCB  separation.  This  means 
that  a  vehicle  which  is  stable  when  properly  trimmed  will  remain  stable 
for  off-trim  conditions. 

2.  Loss  of  stability  occurs  always  in  the  form  of  supercritical  Hopf  bifurca¬ 
tions  with  the  generation  of  stable  limit  cycles.  It  was  found  that  this  is 
mainly  due  to  the  stabilizing  effects  of  the  quadratic  drag  forces. 
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3.  Even  though  the  quadratic  drag  forces  do  not  influence  the  initial  loss  of 
stability,  they  have  a  significant  effect  on  post-loss  of  stability  stabiliza¬ 
tion. 

4.  In  general,  longer  aft  body  sections  seemed  to  increase  the  range  of  linear 
stability  but  influence  adversely  the  resulting  limit  cycles  upon  the  initial 
loss  of  stability. 

It  should  be  emphasized  that  the  occurrence  of  supercritical  Hopf  bifurcations 
is  an  attribute  of  the  open  loop  system  only.  Under  closed  loop  control,  it  is 
possible  to  experience  either  supercritical  or  strongly  subcritical  Hopf  bifurca¬ 
tions,  as  shown  in  [Papoulias  et  al  (1995)].  The  latter  are  particularly  severe 
in  practice  since  self-sustained  vehicle  oscillations  may  be  initiated  prior  to 
loss  of  stability,  depending  on  the  level  of  external  excitation  or  the  initial 
conditions. 
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APPENDIX 


The  following  is  a  list  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN  or  MATLAB.  Complete  print¬ 
outs  of  the  programs  follow  after  the  list. 


•  CRIT_O.M 

MATLAB  program  for  calculating  the  critical  speed  for  (5  =  0. 

•  DSTAB.M 

MATLAB  program  for  calculating  the  degree  of  stability. 

•  HOPF-O.FOR 

FORTRAN  program  for  evaluation  of  hopf  bifurcation  formulas  using  the 
suboff  submarine  model. 
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'/.  Program  crit_0.m 

7o  Evaluation  of  critical  speed  for  delta=0 
clear 

rho  =  1.94; 
g  =  32.2; 

L  =  13.9792; 

ndl  =  0.5*rho*L“2; 

nd2  =  0.5*rho*L*3; 

nd3  =  0.5*rlio*L~4; 

nd4  =  0.5*rho*L“5; 

m  =  1556. 2363/ (g*nd2); 

md  =  1556.2363/g; 

V  =  md/rho ; 

zg  =  0.005; 

while  zg<0.026, 
flagl  =  0; 

for  Fn  =  0.1:0.01:0.35; 
flagl  =  f lagl+1 ; 

flag2  =  0; 

for  Fm  =  0.3:0.01:0.6; 
flag2  =  flag2+l; 

Fb  =  1-Fn-Fm; 

d  =  ((12*V) ./(pi*L*(3*Fm+2*Fn+Fb))) .*0.5; 

r  =  d/2; 

Vn  =  (2/3*pi*r.“2*L.*Fn) ; 

Mn  =  Vn*rho ; 

Vm  =  (pi*r.“2.*Fm*L) ; 

Mm  =  Vm=('rho ; 

Vb  =  (l/3*pi*r.“2*L.*Fb) ; 

Mb  =  Vb*rho ; 

In  =  Mn.*(l/5*(r.~2+(L*Fn) ."2)-(3*L*Fn/8) .*2); 

Im  =  Mm/12.*(3*r.*2+(L*Fm) .*2) ; 

Ib  =  Mb.*(3/5*(r.-2/4+(L=t=Fb)  .~2)-(L*Fb/4)  .*2)  ; 

xcb  =  pi*d.“2.*(2*L=t'Fn.*(L*Fm/2+3=i'L*Fn/8)  . .  . 

-L*Fb.*(L*Fb/4+L*Fm/2))/(12*V) ; 

Lcb  =  L*(Fn+Fm/2)-xcb; 

lyd  =  In+Im+Ib+(Mn. *(Lcb-5*L*Fn/8) . *2)  . . . 

+(Mm. *(Lcb-L*Fm/2-L*Fn) . *2) . . . 
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+  (Mb .  *  (Lcb-L*  (Fn+Fin+Fb/4)  )  .  ~2)  ; 
lym  =  Iyd/nd4; 

%  inputs  A1,A2,A3,A4,A5,A6,A7,A8  for  each  coefficient 
Al= [-0.0641,  0.0277,  -0.0314,  -0.0003,  0.0002,  -0.0002,  -0.0031] 

A2= [-0.1149,  0.0499,  -0.0559,  0.0040,  0.0007,  -0.0007,  -0.0046] 

A3= [-0.0632,  0.0266,  -0.0292,  0.0027,  0.0007,  -0.0007,  -0.0021] 

A4=[  0.0670,-0.0283,  0.0310,  -0.0012,  -0.0008,  0.0008,  0.0031] 

A5=[  0.0732,-0.0301,  0.0316,  -0.0045,  -0.0016,  0.0016,  0.0024] 

A6= [-0.0263, -0.0056,  -0.0091,  0.0006,  -0.0144,  0.0144,  -0.0013] 
A7= [-0.5769, -1.6357,  -0.0880,  -0.1590,  -1.8067,  1.8067,  -0.0808] 


’/,  Hydrodynamic  coefficient  prediction  equation 
Cl  =  8.023e-3; 

for  i  =  1:7, 

HCm(i)  =  Al(i)t'Fn.‘2+A2(i)*Fn.*Fm.  .  . 

+A3(i)*Fm.-2+A4(i)*Fn. .. 
+A5(i)*Fm+A6(i)+A7(i)*(V/L'3-Cl); 

end 

zqdot  =  -6.33e-4; 

HCm(8)  =  zqdot; 

ratio  =  [0.5686,-1. 4357 , -0 . 2658 , 0 . 2675 ,1.1781,-30.5114,0.8149,1.0] 
HC  =  HCm. /ratio; 

zqdot  =  -6.33e-4; 
zwdot  =  HC(5); 
zq  =  HC(3); 

zw  =  HC(1); 

mqdot  =  HC(7); 
mwdot  =  HC(6) ; 
mq  =  HC(4); 

mw  =  HC(2); 

Iratio  =  0.92943; 

ly  =  lym/Iratio; 

cd  =0.015; 

zb  =  0/L; 

xudot  =  -0.05*m; 
xb  =  0/L; 

xg  =  0; 


49 


Gv  =  1  -  mw .  *  (zq+m)  .  /  (zw .  *  (mq-m .  *xg) ) ; 

xgb  =  xg-xb; 
zgb  =  zg-zb; 

for  j  =  1: length (zg) 
for  i  =  l:length(xg) 

theta(i,j)  =  atan(-xgb(i) ./zgb(j) ) ; 

aO  =  (m-zwdot)*(Iy-mqdot)-(mwdot+in*xg(i))*(zqdot+m*xg(i)) ; 
bO  =  (-zwdot*m-m*mw-zq*m)*xg(i)  +  (-m*mq+zwdot*niq-zqdot=t=mw. .  . 

-zq*mwdot-m*inwdot-Iy*zw+mqdot*zw) ; 
cO  =  -in*zw*xg(i)+mq*zw-zq=i'mw-m*mw; 

cl  =  (-m*xg(i)+zwdot*xg(i)+ni*xb-zwdot*xb)*sin(theta(i, j)) . . 

+ (-m*zb-zwdot*zg( j ) +zwdot*zb+m*zg ( j ) ) *cos (thetaCi , j ) ) 
dl  =  (zw*xg(i)-zw*xb)*sin(theta(i, j)) . . . 

+ (zw*zb-zw*zg ( j ) ) *cos (thetaCi , j ) ) ; 

w(i,j)  =  bO*cO/(aO*dl-bO*cl) ; 
uO(i,j)  =  (1556.2363/(ndl*w(i,j)))''.5; 

ucr(flag2,flagl)  =  uO(i,j); 

end 

end 

end 

end 

Fn  =0.1:0.01:0.35; 

Fm  =  0.3:0.01:0.6; 

mesh(Fn,Fm,ucr/8) ,grid 
xlabel(’Fn’) 
ylabel(’Fm’) 
zlabeK  ’ucr’ ) 
hold  on 

zg=zg+0 . 01 ; 

end 
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%  Program  dstab.m 

7o  Matlab  program  for  calculation  the  degree  of  stability 
clear 

clear  global 
rho  =  1.94; 
g  =  32.2; 

L  =  13.9792; 

ndl  =  0 . 5*rho*L'"2 ; 

nd2  =  O.S+rho+L'S; 

nd3  =  0.5*rho*L‘"4; 

nd4  =  O.S^rho’t'L^S; 

m  =  1556. 2363/ (g*nd2); 

md  =  1556.2363/g; 

V  =  md/ rho ; 

fig  =1; 

for  Fn  =  0.1:0. 2:0. 3, 
for  Fm  =  0.4:0. 2:0. 6, 

Fb  =  l-Fn-Fm; 

d  =  ((12*V)  ./(pi*L*(3*Fm+2*Fn+Fb)))  .^'O.b; 

r  =  d/2; 

Vn  =  (2/3*pi*r . "2*L.*Fn) ; 

Mn  =  Vn*rho ; 

Vm  =  (pi*r."2.*Fm*L) ; 

Mm  =  Vm*rho ; 

Vb  =  (l/3*pi*r."2*L.*Fb) ; 

Mb  =  Vb*rho; 

In  =  Mn.  *  (1/5*  (r .  ^2+(L*Fn) .  -'2)-(3*L*Fn/8) .  ^2) ; 

Im  =  Mm/12.*(3*r.^2+(L*Fm) .^2); 

Ib  =  Mb . * (3/5* (r . ^2/4+(L*Fb) . ^2) - (L*Fb/4) . ^2) ; 
xcb  =  pi*d. '"2.*(2*L*Fn.*(L*Fm/2+3*L*Fn/8)  .  .  . 

-L*Fb . * (L*Fb/4+L*Fm/2) ) / (12*V) ; 

Lcb  =  L*(Fn+Fm/2)-xcb; 

lyd  =  In+Im+Ib+(Mn.*(Lcb-5*L*Fn/8) .^2) . . . 

+(Mm. *(Lcb-L*Fm/2-L*Fn) . *2) . . . 
+(Mb.*(Lcb-L*(Fn+Fm+Fb/4)) .^2); 


51 


lym  =  Iyd/nd4; 


7.  inputs  A1,A2 
Al= [-0.0641,  0 
A2= [-0.1149,  0 
A3= [-0.0632,  0 
A4=[  0.0670,-0 
A5=[  0.0732,-0 
A6= [-0.0263,-0 
A7= [-0.5769,-1 


A3,A4,A5,A6,A7 
0277,  -0.0314, 
0499,  -0.0559, 
0266,  -0.0292, 
0283,  0.0310, 
0301,  0.0316, 
0056,  -0.0091, 
6357,  -0.0880, 


A8  for  each 
-0.0003,  0 

0.0040,  0 

0.0027,  0 

-0.0012,  -0 
-0.0045,  -0 
0.0006,  -0 
-0.1590,  -1 


coefficient 
0002,  -0.0002, 
0007,  -0.0007, 
0007,  -0.0007, 
0008,  0.0008, 
0016,  0.0016, 
0144,  0.0144, 

8067,  1.8067, 


-0.0031] 
-0.0046] 
-0.0021] 
0.0031] 
0 . 0024] 
-0.0013] 
-0 . 0808] 


7.  Hydrod3mamic  coefficient  prediction  equation 
Cl  =  8.023e-3; 
for  i=l:7, 

HCm ( i ) =A 1 ( i ) *Fn . ' 2+A2 ( i ) *Fn . *Fm+A3 ( i ) *Fm . *  2+A4 ( i ) *Fn .  . . 
+A5(i)*Fm+A6(i)+A7(i)*(V/L*3-Cl); 


end 

zqdot  =  -6.33e-4; 

HCm(8)  =  zqdot; 

ratio  =  [0.5686,-1.4357,-0.2658, 
HC=HCm. /ratio; 

zqdot  =  -6.33e-4; 
zwdot  =  HC(5) ; 
zq  =  HC(3); 
zw  =  HC(1); 
mqdot  =  HC(7) ; 
mwdot  =  HC(6); 
mq  =  HC(4); 
mw  =  HC(2); 

Iratio  =  0.92943; 
iy  =  lym/Iratio; 

cd  =  0.015; 


.2675,1.1781,-30.5114,0.8149,1.0]  ; 
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zg 

zb 


=  0.015; 

=  0/L; 

m  =  1556, 2363/ (g*nd2); 

xudot  =  “0.05*m; 

xb  =  0/L; 

xg  =  linspace (~0 .01,0.01,41); 

uo  =  8*linspace (0.2, .6,41) ; 

Gv  =  1  -  inw.*(zq+m)  ./(zw.*(mq-m.*xg)) ; 

w  =  1556.2363./(ndl.*uo.’'2)  ; 

b  =  w; 

xgb  =  xg-xb ; 

zgb  =  zg-zb ; 

theta  =  at an (“Xgb. /zgb ) ; 
for  j  =  l:length(uo) 
for  i  =  1: length (xg) 

A  =  [-2*cd  00  0; 

0  zw  (zq+m)  0 ; 

0  mw  (mq-in*xg(i))  (xgb(i)*sin(theta(i)) . .  . 

-zgb*cos (theta(i) ) ) *b( j ) 
0  0  1  0]  ; 

B  =  [m-xudot  0  ni*zg  0; 

0  m-zwdot  -(m*xg(i)+zqdot)  0; 

in*zg  -(mwdot+ni*xg(i))  iy-mqdot  0; 

0  0  0  1]  ; 

evalsl  =  eig(A,B);  %  no  surge  coupling 
degstabKi,  j)  =  mcLX (real (evalsl) )  ; 
end 
end 

figure (fig) 

mesh(uo/8,xg,degstabl) ,grid 
xlabeK^uoO 
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ylabel(’xg’) 

zlabeK’ degree  of  stability’) 
fig=fig+l; 
end 
end 
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C  PROGRAM  HOPF.FOR 

C  EVALUATION  OF  HOPF  BIFURCATION  FORMULAS 

C  USING  THE  SUBOFF  SUBMARINE  MODEL 

C  Cd=0.5,  ZG=0.015  ,Fn=0.24,  Fm=0.52 


1 

2 

3 

4 

5 

6 
7 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  L, IY,MASS,MQDOT,MWDOT,NDl , 
MQ,MW,K1,K2, 

BETA, GAMA, 

E0,E1,E2,E3,E4, 

DW1,DW2,DW3,DV4, 

DQl,Dq2,DQ3,DQ4, 

MASSM , MASSN , MASSB , IB , IM , IN , 

RHO , CD , RADI , VOLM , VOLN , VOLB , LCB , XCB 


C 


DOUBLE  PRECISION 

1 

2 

3 

4 

5 

6 
7 

DOUBLE  PRECISION 


M11,M12,M13,M21,M22,M23, 

M31,M32,M33, 

N11,N12,N13,N21,N22,N23, 

N31,N32,N33, 

L2 1 , L22 , L23 , L24 , L31 , L32 , L33 , L34 , 
L25 , L26 , L27 , L35 , L36 , L37 , 

L2 1 A , L22A , L23A , L24A , L3 lA , 

L32A,L33A,L34A 

LN,LM,LB,FM,FN,FB,KK 


DIMENSION  A(3,3) ,T(3,3) ,TINV(3,3) ,FV1(3) , IVl (3) ,YYY(3,3) 
DIMENSION  WR(3) ,WI(3) ,TSAVE(3,3) ,TLUD(3,3) ,IVLUD(3) 
DIMENSION  ASAVE(3,3) ,A2(3,3) ,XL(55) ,BR(55) 

DIMENSION  VEC0(55),VEC1(55) ,VEC2(55) ,VEC3(55) ,VEC4(55) 
DIMENSION  HCA1(7) ,HCA2(7) ,HCA3(7) ,HCA4(7) ,HCA5(7) 
DIMENSION  HCA6(7) ,HCA7(7) ,HC(8) ,RATI0(7) , SVLUD(3) 


C 

OPEN  (20 , FILE= ’ DATA.O . DAT  ^  STATUS= ’ NEW  O 
C 
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WEIGHT=  1556.2363 

L  =  13.9792 

RHO  =  1.94 

DO  8886  GDI  =  0.40,0.60,0.10 

CD  =  0.5^'CD1*RH0 

G  =  32.2 

XB  =  0.0 

DO  8887  KK  =  0.0050,0.0250,0.0050 

ZG  =  KK*L 

ZB  =  0.0 

MASS  =  WEIGHT/G 

BOY  =  WEIGHT 

VOLUME=  MASS/RHO 

DO  8888  FN=0. 10,0.32,0. 10 

DO  8889  FM=0. 40, 0.62, 0.10 


C  WRITE  (20,*)  ’CD  =’ ,CD 
C  WRITE  (20,*)  ’ZG  =’ ,KK 
C  WRITE  (20,*)  ’FN  =’ ,FN 
C  WRITE  (20,*)  ’FM  =’,FM 


FB  =  1.0-FN-FM 

LN  =  L*FN 

LM  =  L*FM 

LB  =  L*FB 

DIAM  =  SQRT((12.*V0LUME) 

&  /(3.14159*L*(3.*FM+2.*FN+FB))) 

WRITE(*,4001)  DIAM 
RADI  =  DIAM/2. 

VOLN  =  (2./3.*3.14159*RADI**2.*L*FN) 

MASSN  =  VOLN*RHO 

VOLM  =  (3.14159*RADI**2.*FM*L) 

MASSM  =  VOLM*RHO 

VOLB  =  (l./3.*3.14159*RADI**2.*L*FB) 

MASSB  =  VOLB*RHO 

IN  =  MASSN* (l./5.*(RADI**2+(L*FN)**2) 
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&  -(3*L*FN/8)**2) 

IM  =  MASSM/12.*(3.*RADI**2.  +  (L*FM)*=t=2) 

IB  =  MASSB*(3./5.*(RADI**2/4. 

&  +a*FB)**2)-a*FB/4.)**2) 

XCB  =  (3.14159*DIAM**2*(2.*L*FN 

&  *(L*FM/2.+3.*L*FN/8.) 

&  -L*FB* (L*FB/4 . +L*FM/2 .)))/( 12 . *VOLUME) 

WRITE (*,4001)  XCB 
LCB  =  L*(FN+FM/2.)-XCB 

WRITE(*,4001)  LCB 

lY  =  IN+IM+IB+MASSN*(LCB-5*L*FN/8)**2 

&  +MASSM* (LCB-L*FM/2-L*FN) **2 

lY  =  IY+MASSB*(LCB-L*(FN+FM+FB/4))**2 


C  Inputs  A1,A2,A3,A4,A5,A6,A7,A8  for  each  coefficient 


HCAl(l) 

= 

-0.0641 

HCA1(2) 

= 

0.0277 

HCA1(3) 

= 

-0.0314 

HCA1(4) 

= 

-0.0003 

HCA1(5) 

= 

0.0002 

HCA1(6) 

= 

-0 . 0002 

HCAl (7) 

= 

-0.0031 

HCA2(1) 

= 

-0.1149 

HCA2(2) 

= 

0 . 0499 

HCA2(3) 

= 

-0.0559 

HCA2(4) 

= 

0.0040 

HCA2(5) 

= 

0.0007 

HCA2(6) 

= 

-0.0007 

HCA2(7) 

= 

-0.0046 

HCA3(1) 

= 

-0.0632 

HCA3(2) 

= 

0.0266 

HCA3(3) 

= 

-0.0292 

HCA3(4) 

= 

0.0027 

HCA3(5) 

= 

0.0007 

HCA3(6) 

= 

-0.0007 
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HCA3(7)  =  -0 

HCA4(1)  =  0 

HCA4(2)  =  -0 

HCA4(3)  =  0 

HCA4(4)  =  -0 

HCA4(5)  =  -0 

HCA4(6)  =  0 

HCA4(7)  =  0 

HCA5(1)  =  0 

HCA5(2)  =  -0 

HCA5(3)  =  0 

HCA5(4)  =  -0 

HCA5(5)  =  -0 

HCA5(6)  =  0 

HCA5(7)  =  0 

HCA6(1)  =  -0 

HCA6(2)  =  -0 

HCA6(3)  =  -0 

HCA6(4)  =  0 

HCA6(5)  =  -0 

HCA6(6)  =  0, 

HCA6(7)  =  -0, 

HCA7(1)  =  -0, 

HCA7(2)  =  -1, 

HCA7(3)  =  -0, 

HCA7(4)  =  -0. 

HCA7(5)  =  -1. 

HCA7(6)  =  1. 

HCA7(7)  =  -0. 

C  Hydrodynamic 

C 

Cl 

RATIO (1)  = 

RATIO (2)  =  - 

RATIO (3)  = 


.0021 

.0670 

.0283 

.0310 

.0012 

.0008 

.0008 

.0031 

.0732 

.0301 

.0316 

.0045 

.0016 

.0016 

.0024 

.0263 

.0056 

.0091 

.0006 

.0144 

0144 

0013 

5796 

6357 

0880 

1590 

8067 

8067 

0808 


coefficient  prediction  equation 


=  8.023E-03 

0.5686 
1.4357 
-0.2658 
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RATIO (4)  =  0.2675 

RATIO (5)  =  1.1781 

RATIO (6)  =  -30.5114 

RATIO (7)  =  0.8149 

DO  5000  1=1,7 

HC ( I ) = (HCAl ( I ) *FN**2+HCA2 ( I ) *FN*FM 
&  +HCA3(I)*FM**2+HCA4(I)*FN 

&  +HCA5(I)*FM+HCA6(I) 

&  +HCA7(I) * (VOLUME/ (L*L*L) -Cl) ) /RATIO (I) 


5000  CONTINUE 

HC(8)  =  -6.33E-04 

ZQDOT  =  -6 . 33E-04*0 . 5*RH0*L**4 

HC(8)  =  ZQDOT 

ZWDOT  =  HC ( 5) *0 . 5*RH0*L*  *3 

ZQ  =  HC(3)*0.5*RH0*L**3 

ZW  =  HC(1)*0.5*RH0*L**2 

MQDOT  =  HC(7)*0.5*RH0=t=L=t'*5 

MWDOT  =  HC(6)*0.5*RH0*L=i'*4 

MQ  =  HC(4)*0.5*RHO»L**4 

MW  =  HC(2)*0.5*RH0=t=L**3 

RATIO 1  =  0.92943 

lY  =  lY/RATIOl 

WRITE (*,4001)  lY 

NDl  =  0.5*RH0*L**2 

ZGB  =  ZG-ZB 

C 

C  DEFINE  THE  LENGTH  X  AND  BREADTH  B  TERMS  FOR  THE  INTEGRATION 
C 


DO  333  1=0,21 
XL(I+1)=  I*LB/21.0 
BR(I+1) =DIAM*XL(I+1) /LB 
333  CONTINUE 
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DO  334  1=1,2 
XL(22+I)=  LB+I*LM/2.0 
BR(22+I)=DIAM 

334  CONTINUE 

DO  335  1=1,30 
c  WRITE(*,*)  I 

XL (1+24) =  XL  Cl+23) +1 . /4 . * (L-XL (1+23) ) 

IF  (((XL(I+24)-LB-LM)**2/(LN**2)).GT.1.0)  THEN 
BR(I+24)=0.0 
ELSE 

BR(I+24) =DIAM*SQRT( 1 . 0- ( (XL(I+24) -LB-LM) **2/ (LN**2) ) ) 
ENDIF 

335  CONTINUE 

XL (55)  =  L 
BR(55)  =  0 
DO  102  N  =  1,55 

XL(N)  =  XL(N)-L+LCB 
102  CONTINUE 

WRITE (20, 7001)  XL 
WRITE(20,7001)  BR 
C 

DO  104  K  =  1,55 
VECO(K)=BR(K) 

VEC1(K)=XL(K)*BR(K) 

VEC2(K)=XL(K)*XL(K)*BR(K) 

VEC3 (K) =XL (K) *XL (K) *XL (K) *BR (K) 

VEC4 (K) =XL (K) *XL (K) *XL (K) *XL (K) *BR (K) 

104  CONTINUE 

CALL  TRAP(55,VEC0,XL,E0) 

CALL  TRAP(55,VEC1,XL,E1) 

CALL  TRAP(55,VEC2,XL,E2) 

CALL  TRAP(55,VEC3,XL,E3) 

CALL  TRAP(55,VEC4,XL,E4) 

EPSILON  =  0.001 
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XGMIN=-0.01 

XGMAX=+0.01 

IXG=80 

XGMIN=XGMIN*L 

XGMAX=XGMAX*L 


DO  1  IT=1,IXG 

C  WRITE  (*,3001)  IT.IXG 

XG=XGMIN+ (XGMAX-XGMIN) * ( IT- 1 ) / ( IXG- 1 ) 

XGB=XG-XB 

DV= (MASS-ZWDOT) * (lY-MQDOT) 

&  -(MASS*XG+ZQDOT)*(MASS*XG+MWDOT) 

CD6=CD/ (3 . DO*EPSILON*DV) 

DW1=CD6* ( ( lY-MQDOT) * (-E0) + (MASS*XG+ZQDOT) *E1) 

DW2=CT»6* ( (lY-MQDOT) * (3*E1) - (MASS*XG+ZQDOT) *3*E2) 
DW3=CD6* ( ( lY-MQDDT) * (-3*E2) + (MASS*XG+ZQDOT) *3*E3) 
DW4=CD6* ( (lY-MQDOT) *(E3) - (MASS*XG+ZQDOT) *E4) 

DQ1=CD6* ( (MASS-ZWDOT) * (El) + (MASS*XG+MWDOT) * (-E0) ) 
DQ2=CD6* ( (MASS-ZWDOT) * (-3*E2) +(MASS*XG+MWDOT) * (3*E1) ) 
DQ3=CD6* ( (MASS-ZWDOT) * (3*E3) + (MASS*XG+MWDOT) * (-3*E2) ) 
DQ4=CD6* ( (MASS-ZWDOT) * (-E4) + (MASS*XG+MWDOT) * (E3) ) 
THETAO=ATAN(-XGB/ZGB) 

AAO= (MASS-ZWDOT) * ( I Y-MQDOT) 

&  - (MWDOT+MASS*XG) * (ZQDOT+MASS*XG) 

BBO= (-ZWDOT*MASS-MASS*MW-ZQ*MASS) *XG 
&  +(-MASS*MQ+ZWDOT*MQ-ZClDOT*MW 

&  -ZQ*MWDOT-MASS*MWDOT-IY*ZW+MQDOT*ZW) 

CCO=-MASS*ZW*XG+MQ*ZW-ZC)*MW-MASS*MW 
CC1=(-MASS*XG+ZWDOT*XG+MASS*XB-ZWDOT*XB)*SIN(THETAO) 
&  + (-MASS*ZB-ZWDOT*ZG+ZWDOT*ZB+MASS*ZG) *COS (THETAO) 

DD1=(ZW*XG-ZW*XB)*SIN(THETA0)+(ZW*ZB-ZW*ZG) 

&  *COS (THETAO) 


C  After  applying  AD=BC  (  Routh  Criterion  ) ,  we  manage  to  calculate 
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C  the  critical  speed  UO. 


WEI=BBO*CCO/ (AAO*DDl-BBO*CCl ) 

UO=DSQRT ( 1556 . 2363/WEI) 

U=U0 

WRITE  (*,*)  U/8.0,XG/L 
C 

C  DETERMINE  [A]  AND  [B]  COEFFICIENTS 

C 

A11DV=  (lY-MQDOT)  *ZW+  (MASS*XG+ZC1D0T)  *MW 
A12DV= (lY-MQDOT) * (MASS+ZQ) +(MASS*XG+ZQDOT) * (MQ-MASS*XG) 
A13DV=- (MASS*XG+ZQDOT) ^WEIGHT 
A21DV= (MASS-ZWDOT) ♦MW+ (MASS*XG+MWDOT) *ZW 
A22DV= (MASS-ZWDOT) * (MQ-MASS*XG) + (MASS*XG+MWDOT) * (MASS+ZQ) 
A23DV=- (MASS-ZWDOT) ♦WEIGHT 
C 

A11=A11DV/DV 

A12=A12DV/DV 

A13=A13DV/DV 

A21=A21DV/DV 

A22=A22DV/DV 

A23=A23DV/DV 

C 

C11DV=(IY-MQD0T)^MASS^ZG 
C12DV=- (MASS^XG+ZQDOT) ♦MASS^ZG 
C21DV=- (MASS-ZWDOT) ♦MASS^ZG 
C22DV= (MASS^XG+MWDOT) ♦MASS^ZG 
C 

C11=C11DV/DV 

C12=C12DV/DV 

C21=C21DV/DV 

C22=C22DV/DV 

C  EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 

C 

Kl=- (XGB^SIN (THETAO) -ZGB^COS (THETAO) ) 
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K2=- ( 1 . /6 . ) * (ZGB*COS (THETAO) -XGB*SIN (THETAO) ) 


A(1,1)=0.0 
A(l,2)=0.0 
A(l,3)=1.0 
A(2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1 
A(3,2)=A21*U 
A(3,3)=A22*U 
DO  11  1=1,3 
DO  12  J=l,3 

ASAVE(I,J)=A(I,J) 

12  CONTINUE 
11  CONTINUE 

CALL  RG(3,3,A,WR,WI,1,YYY,IV1,FV1,IERR) 
CALL  DSOMEG ( lEV , WR , WI , OMEGA , CHECK) 

C  WRITE  (*,*)  lEV 

C  WRITE  (*,*)  (WR(IWR),IWR=1,3) 

C  WRITE  (*,*)  (WI(IWI),IWI=1,3) 

OMEGAO=OMEGA 
DO  5  1=1,3 

T(I,1)=  YYY(I,IEV) 
T(I,2)=-YYY(I,IEV+1) 

5  CONTINUE 

IF  (lEV.EQ.l)  GO  TO  13 
IF  (IEV.EQ.2)  GO  TO  14 
STOP  3004 
14  DO  6  1=1,3 

T(I,3)=YYY(I,1) 

6  CONTINUE 
GO  TO  17 

13  DO  16  1=1,3 

T(I,3)=YYY(I,3) 

16  CONTINUE 


63 


CONTINUE 


NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 

INORM=l 

IF  (INORM.NE.O)  CALL  NORMAL(T) 

INVERT  TRANSFORMATION  MATRIX 


DO  2  1=1,3 
DO  3  J=l,3 
TINV(I,J)=0.0 
TSAVECI, J)=T(I, J) 

CONTINUE 

CONTINUE 

CALL  DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO  4  1=1,3 

IF  (IVLUD(I) .EQ.O)  STOP  3003 
CONTINUE 

CALL  DILU(3,3,TLUD,IVLUD,SVLUD) 

DO  8  1=1,3 
DO  9  J=l,3 

TINV(I,J)=TLUD(I,J) 

CONTINUE 

CONTINUE 

CHECK  Inv(T)*A*T 


IMULT=1 

IF  (IMULT.EQ.l) 
IF  (IMULT.EQ.O) 
P=A2(3,3) 

PEIG=P 

VfRITE  (*,4001) 
WRITE  (*,4001) 
WRITE  (*,4001) 


CALL  MULT (TINV, ASA VE,T,A2) 
STOP 


(A2(1,JA2),JA2=1,3) 
(A2(2,JA2),JA2=1,3) 
(A2(3,JA2) ,JA2=1,3) 


PAUSE 


C 
C 

C  DEFINITION  OF  Nij 

C 

N11=TINV(1,1) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N13=TINV(1,3) 

N23=TINV(2,3) 

N33=TINV(3,3) 

C 

C  DEFINITION  OF  Mij 

C 

M11=T(1,1) 

M21=T(2,1) 

M31=T(3,1) 

M12=T(1,2) 

M22=T(2,2) 

M32=T(3,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 

C 

C  DEFINITION  OF  Lij 

C 

L25=C11*M31*M31+C12*M21*M31 
L26=2*C11*M31*M32+C12*(M21*M32+M22*M31) 
L27=C1  l*M32=i'M32+C12*M22=t=M33 
L35=C22*M31*M31+C21*M21*M31 
L36=2*C22*M31*M32+C21*(M21*M32+M22*M31) 
L37=C22=t=M32*M32+C21*M33=t=M22 
C 

C  DEFINITION  OF  ALFA,  BETA,  GAMA 
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D1  =N32*L25  +  N33*L35 
D2  =N32*L26  +  N33*L36 
D3  =N32*L27  +  N33*L37 

Dll=-P 

D12=0MEGA0 

D21=-2*0MEGA0 

D22=-P 

D23=2*0MEGA0 

D32=-0MEGA0 

D33=-P 


BETA=(D2-D21*D1/D11-D23*D3/D33) 

&  / (022-021*012/011-023*032/033) 

ALFA= (D1-012*BETA) /Oil 
GAMA=(03-032*BETA) /033 

L21A=2*C11*ALFA*M31*M33+C12*ALFA 

*(M21*M33+M23*M31) 

L22A=2*C11*ALFA*M32*M33  +  2*C11*BETA*M31*M33 
+  C12*ALFA*(M22*M33+M32*M23) 

+  C12*BETA*(M21*M33+M23*M31) 

L23A=2*C11*GAMA*M31*M33+2*C11*BETA*M32*M33 
+  C12*GAMA*(M21*M33+M23*M31) 

+  C12*BETA*(M22*M33+M23*M32) 

L24A=2*C11*GAMA*M32*M33+C12*GAMA 
&  * (M22*M33+M23*M32) 

L31A=2*C22*ALFA*M31*M33+C21*ALFA 
&  * (M2 1 *M33+M23*M3 1 ) 


L32A=2*C22*ALFA*M32*M33+2*C22*BETA*M31*M33 


+  C21*ALFA*(M22*M33+M32*M23) 

+  C21*BETA*(M21*M33+M23*M31) 

L33A=2*C22*GAMA*M31*M33+2*C22*BETA*M32*M33 
+  C21*GAMA*(M21*M33+M23*M31) 

+  C21*BETA*(M22*M33+M23*M32) 

L34A=2*C22*GAMA*M32  *M33+C2 1 *GAMA 
*(M22*M33+M23*M32) 

L21=L21A+A13*K2*M11**3+DW1*M21**3 
+DW2*M31*M21**2 
+DW3*M21*M31**2+DW4*M31**3 
L22=L22A+3*A13*K2*M12*M11**2+3*DW1*M22*M21**2 
+DW2*  (2*M21*M22*M31+M32*M21*=t=2) 
+DW3* (2*M21*M31*M32+M22*M31**2) 
+3*DW4*M32*M31**2 

L23=L23A+3*A13*K2*Mll*M12*=i'2+3*DWl*M21*M22**2 
+DW2* (M31*M22**2+2*M21*M22*M32) 
+  DW3*(M21*M32**2+2*M22*M31*M32) 

+  3*DW4*M31*M32**2 

L24=L24A+A13*K2*M12**3+DW1*M22**3 
+DW2*M32*M22*=t:2 
+DW3*M22*M32**2+DW4*M32**3 
L31=L31A+A23*K2*M11**3+DQ1*M21**3 
+DQ2*M31*M21**2 
+DQ3*M21*M31**2+DQ4*M31**3 
L32=L32A+3=t'A23*K2*M12*Mll*=K2+3*DQl*M22*M21**2 
+DQ2* (2*M21*M22*M31+M32*M21**2) 
+DQ3* (2*M21*M31*M32+M22*M31**2) 
+3*DQ4*M32*M31**2 

L33=L33A+3*A23*K2*M11*M12**2+3*DQ1*M21*M22**2 
+  DQ2*(M31*M22**2+2*M21*M22*M32) 

+  DQ3*(M21*M32*=t'2+2*M22*M31*M32) 

+  3*DQ4*M31*M32**2 

L34=L34A+A23*K2*M12*=t=3+DQl*M22i'*3 


&  +DQ2*M32*M22**2 

&  +DQ3*M22*M32**2+DQ4*M32**3 

C 

R11=N12*L21+N13*L31 

R12=N12*L22+N13*L32 

R13=N12*L23+N13*L33 

R14=N12*L24+N13*L34 

R21=N22*L21+N23*L31 

R22=N22*L22+N23*L32 

R23=N22*L23+N23*L33 

R24=N22*L24+N23*L34 

C 

C  EVALUATE  DALPHA  AND  DOMEGA 

C 

UINC=0.001 
UR  =U+UINC 
UL  =U-UINC 
U  =UR 
C 

A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22=t=U 

C 

CALL  RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREC) 

C 

U=UL 

C 
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A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 

C 

CALL  RG(3,3,A,WR,WI,0,YYy,IVl,FVl,IERR) 

CALL  DSTABL (DECS, WR,WI, FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 

C 

DALPHA= (ALPHR-ALPHL) / (UR-UL) 

DOMEGA= (OMEGR-OMEGL) / (UR-UL) 

C 

C  EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 

C 

C0EF1=3 . 0*R1 1+R13+R22+3 . 0*R24 
C0EF2=3 . 0*R21+R23-R12-3 . 0*R14 
AMU2  =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*C0EF1 

C  TAU2  =-(COEF2-DOMEGA*COEFl/DALPHA)/(8.0*OMEGAO) 

C  PER  =2.0*3. 1415927/OMEGAO 

C  PER  =PER*U/L 

C  WRITE  (20,2001)  XCB 

WRITE  (20,2001)  XG/L,C0EF1 
1  CONTINUE 
8889  CONTINUE 
8888  CONTINUE 
8887  CONTINUE 
8886  CONTINUE 
STOP 

1001  FORMAT  (’  ENTER  NUMBER  OF  DATA  LINES’) 
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1002  FORMAT 

1003  FORMAT 

1004  FORMAT 
2001  FORMAT 
4001  FORMAT 
7001  FORMAT 
3001  FORMAT 

END 


(»  ENTER  UO,  ZG,  AND  DSAT’) 

(’  ENTER  BOW  PLANE  TO  STERN  PLANE  RATIO’) 
(’  ENTER  ZG’) 

(2E14.5) 

(1F15.5) 

(6F15.5) 

(215) 
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